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Q ■ Abstract. We study the occurrence of physical collisions between 

stars in young and compact star cluster. The calculations are per- 
formed on the GRAPE-4 with the starlab software environment 
which include the dynamical evolution and the nuclear evolution of 
all stars and binaries. The selection of the initial conditions is based 
t^. ; on existing and well observed star clusters, such as R136 in the 30 

Doradus region in the Large Magellanic Cloud and the Arches and 
Quintuplet star clusters in the vicinity of the Galactic center. Colli- 
sions between stars occurred rather frequently in our models. At any 
time a single star dominates the collision history of the system. The 
collision rate of this runaway merger scales with the initial relaxation 
& ■ time of the cluster and is independent on other cluster parameters, 

^! such as the initial mass function or the initial density profile of the 

cluster. Subsequent encounters result in a steady grow in mass of 
the coagulating star, until it escapes or explodes in a supernova. The 
c3 \ collision rate in these models is about 2.2 x 10~ 4 collisions per star 

per Myr for a cluster with an initial relaxation time of 1 Myr. 
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1. Introduction 

The central regions of young star clusters, globular clusters and galaxies 
have a high stellar density and close encounters between stars are thought 
to occur on a regular basis. In some cases stars may even collide and 
form new objects. However, proof for the occurrence of such collisions 
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has never been found explicitly in observations or in theoretical models. 
Direct evidence for collisions could be provided by the observation of a blue 
straggler with a mass exceeding three times the mass of the clusters' turn 
off, by discovery of a blue straggler in an eccentric orbit around another 
main-sequence star or by catching a collision in the act. 

There are quite a few indications that collisions play a major role in the 
evolution of star clusters and their constituents. For example: the presence 
of blue stragglers and millisecond pulsars in globular clusters indicates that 
stellar collisions may play a role, even though alternative -non-collisional- 
scenarios have also been proposed (Leonard 1989). 

Collision rates in star clusters have always been estimated using cross 
section arguments. Pioneering work in this field has been performed by 
Hills (1975) and Lightman & Shapiro (1977). Their general conclusions 
are: 

1) In order to get a runaway growth of a single star one requires a star 
cluster which contains typically more than ~ 10 7 stars. 

2) A massive star formed out of multiple collisions will always be accom- 
panied by many objects which experienced only one or two collisions. 

3) Stellar collisions are so rare that they do not affect the evolution of 
the star cluster. 

We show that these results are based on insufficient detail in the mod- 
els and that the contrary has to be concluded: 

ad 1) Runaway growth is a natural phenomena which occurs even in rather 
small cluster. 

ad 2) The massive product of many collisions will therefore be well sepa- 
rated in mass from the other objects. 

ad 3) Even in star clusters where physical collisions are relatively rare the 
evolution of the star cluster is altered significantly. 

We will address these issues by the direct integration of the orbits of all 
the stars in such clusters, while accounting for the tidal field of the Galaxy 
and the internal evolution of the stars and binaries. Our model calculations 
are performed using the starlab software environment with up to 32k stars. 
The calculations are sped up with the special purpose computer GRAPE- 
4. The initial conditions are taken to represent observed compact star 
clusters such as R136 (the compact star cluster in the 30 Doradus region), 
the Arches cluster and the Quintuplet system. 
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1.1. The model 



To model a dense star cluster we use a hydride computer program that 
consists of two independent parts. An A^-body model integrates the equa- 
tions of motion of all stars and at the same time the other part computes 
the evolution of the stars. The back coupling between the stellar evolution 
and the stellar dynamics is taken into account in a self consistent -object 
oriented- fashion. 

All computations are performed using the starlab toolset (Portegies 
Zwart et al. 2000b, see also |http : / /www . manybody . o"rg| ). Starlab consists 
of the A^-body integrator kira and the stellar and binary evolution program 
SeBa. The calculations are accelerated using the special purpose computer 
GRAPE-4 (Makino et al. 1997). 

2. Initial conditions 

Selection for the initial conditions in our calculations were guided by the 
well studied Arches and the Quintuplet clusters and the central star cluster 
R136 in the 30 Doradus region of the Large Machelanic cloud. A compi- 
lation of the characteristics of these clusters is presented in Tab. 1. The 
Arches and Quintuplet cluster are very close to the Galactic center where 
R136 is completely isolated from the tidal perturbation of the Galaxy. 

Table 1. Observed parameters for R136, the Arches and the 
Quintuplet systems. Columns list cluster name, reference, age, 
mass, projected distance to the Galactic center, tidal radius 
(jtide), and half mass radius (Vh m )- The final column presents 
an estimate of the density within the half mass radius. 



Name 


ref 


Age 


M 


r GC 


r tide 


^hm 


log Pcorc 






[Myr] 


[1O 3 M ] 




- [pc] - 




[M /pc 3 ] 


R136 


a 


2-4 


21-79 


50k 


> 20 


~ 0.5 


4.6-5.2 


Arches 


b 


1-2 


12-50 


30 


1 


0.2 


5.6-6.2 


Quintuplet 


c 


3-5 


10-16 


50 


1 


~ 0.5 


4.3-4.5 



References: a) Brandl et al. (1996); Campbell et al. (1992); Massey & 
Hunter (1998). b) Figer et al. (1999); c) Glass, Catchpole & Whitelock 
(1987); Figer, Mclean & Morris (1999). 



We performed a total of 44 A^-body calculations over a wide range 
of initial conditions. The number of stars was varied from lk (1024) to 
32k. Initial density profiles and velocity dispersion for the models are 
taken from Heggie-Ramamani models (Heggie & Ramamani 1995) with 
Wo ranging from 1 to 7. At birth, the clusters are assumed to perfectly fill 
the zero velocity surface in the tidal field of the Galaxy. In most cases we 
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selected the initial mass function between 0.1 M and 100 M suggested 
for the Solar neighborhood by Scalo (1986), but several calculations are 
performed with power law initial mass functions with an exponent of - 
2 and -2.35 (Salpeter). The specification of the remaining parameters, 
such as the strength of the Galactic tidal field, the virial radius and the 
tidal radius, are discussed in a forthcoming paper. Here we only specify 
the initial relaxation time at the half mass radius of our models as this 
happens to be the fundamental parameter with which our results can be 
understood (see sect. §0). An overview of the initial conditions for the 
computed models is summarized in Tab.§. A detailed discussion about 
several of the presented calculations by Portegies Zwart et al. (2001). 

The A^-body system is fully determined once the tidal field, the initial 
mass function, the number of stars and the relaxation time of the cluster 
are selected. The total mass of the stellar system determines the unit of 
mass in the A^-body system, the tidal radius r t id e sets the distance unit 
and the velocity dispersion together with the size of the stellar system sets 
the time scale (see Heggie & Mathieu 1986). The evolution of the cluster is 
subsequently followed using the direct A^-body integration including stellar 
and binary evolution and the tidal field of the Galaxy (see Portegies Zwart 
et al. 2000b). For economic reasons not all stars are kept in the A^-body 
system, but stars are removed when they are 3r t id e from the center of the 
star cluster. 

The evolution of our star clusters is driven by two-body relaxation, 
by stellar mass loss and by the external tidal field of the parent Galaxy. 
The initial relaxation time, however, appears to be the most fundamen- 
tal parameter for the range of parameters we study here. The two-body 
relaxation time: 



Here A" is the number of stars and 7 is a scaling factor, introduced to 
model the effects of the cut-off in the long range Coulomb logarithm (see 
Giertz & Heggie 1994; 1996). Here t crss is the half-mass crossing time of 
the cluster is 



Here rh m is its half mass radius and M is the mass of the cluster. 
2.1. Scaling the collision cross section 

The rate at which stars in a cluster experience collisions can be estimated 
via 



Here n c is the number density of the stars in the core, a is the collision 
cross section (for approach within some distance d), and v is the velocity 
dispersion of the cluster stars. 





n co n oc n c av. 



(3) 
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Table 2. Overview of the calculations. The first column gives 
the name of the model as used in previous publications (names 
RxWx and KMLx are from Portegies Zwart et al. 2001 [see also 
Portegies Zwart et al. 2000a]; and the other models are described 
in detail by Portegies Zwart et al 1999). The following columns 
give the number of stars (in units of 1024), the initial mass func- 
tion (Scalo 1986, or a power law with slope as indicated), the ini- 
tial King Wq, the initial relaxation time (in Myr) and the number 
of runs performed with these initial conditions (W run ). The last 
three columns give the average number of collisions in these calcu- 
lations, the moment the last collision occurred and the collision 
rate per Myr per star (see Eq. [| The models indicated with * 
have been computed without a Galactic tidal field (see Portegies 
Zwart 1999). 



model 


(N) 


IMF 


(Wo) 


(w) 


A/run 


(A^coll) 


(^last) 


fcoll 


R34W7 


12k 


Scalo 


7 


0.4 


2 


16. 


10.4 


-3.70 


KML112 


4k 


-2 


7 


0.5 


2 


4.0 


1.9 


-3.00 


KML101 


4k 


-2 


4 


1.4 


2 


2.0 


1.0 


-3.60 


KML142 


6k 


-2.35 


4 


1.9 


1 


1.0 


2.2 


-4.13 


KML111 


4k 


-2 


1 


2.3 


2 


0.5 


6.7 


-4.95 


R90W7 


12k 


Scalo 


7 


2.8 


1 


13. 


10.0 


-4.44 


R34W4 


12k 


Scalo 


4 


3.2 


3 


6.3 


30.0 


-4.26 


KML144 


14k 


-2.35 


4 


3.9 


1 


2.0 


2.4 


-4.18 


R150W7 


12k 


Scalo 


7 


4.5 


2 


10. 


21.3 


-4.52 


6k6X5* 


6k 


Scalo 


6 


5.0 


1 


21. 


47.9 


-3.47 


R34W1 


12k 


Scalo 


1 


5.5 


3 


4.7 


29.1 


-4.28 


R90W4 


12k 


Scalo 


4 


8.1 


5 


5.8 


10.0 


-4.47 


Nk6X10* 


9k 


Scalo 


6 


10.0 


8 


10. 


18.0 


-3.94 


R150W4 


12k 


Scalo 


4 


13.0 


4 


8.5 


7.3 


-4.70 


R90W1 


12k 


Scalo 


1 


14.6 


1 


7.0 


9.8 


-4.57 


6k6X20* 


6k 


Scalo 


6 


20.0 


2 


4.0 


95.4 


-4.19 


R150W1 


12k 


Scalo 


1 


23.6 


2 


3.0 


2.1 


-4.87 


R300W4 


12k 


Scalo 


4 


55.6 


1 


1.0 


10.0 


-5.77 


R34W1 


32k 


Scalo 


1 


58.1 


1 


4.0 


35.8 


-5.89 
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The central density and the collision cross section are given by the 
following proportionalities: 

N 

n c oc — , 

' c 

a oc d 2 + 4- (4) 

v 

We neglect the d 2 term in the cross section. In the case of constant crossing 
time, Nr 3 = constant, we may write v oc N l l 3 . The number of collisions 
is then computed by substitution of Eq. (^) into Eq. ([5]) to obtain n co \\ ~ 
(IN" 1 / 3 resulting in a collision rate per star per relaxation time of: 

dN 2 ' 3 

T^coll ^ — • (5) 

''tlx. 

From a theoretical perspective one expects therefore that the collision 
rate scales with the number of stars and inversely with the initial relaxation 
time of the stellar system. The normalization constant, however, is hard 
to estimate from first principals. 



3. Results 

Shortly after the start of our calculations, each model experiences core 
collapse. This event is triggered by the more massive stars in the cluster, 
which sink into the cluster core in a fraction of the initial relaxation time; 
~ t r i x m/(m). Once in the cluster core these massive stars dominate the 
dynamics by forming binaries with other stars. This phase of cluster evo- 
lution continues until the massive star is either kicked out of the cluster 
by a strong encounter or the star explodes in a supernova (see Portegies 
Zwart et al 1999). The period over which the massive star can dominate 
the cluster dynamics depends on the initial conditions of the star cluster. 

All our models finally dissolve in the tidal field of the Galaxy. (Except 
the isolated models which are terminated after one initial relaxation time 
lapsed.) 

Once a massive star settles in the cluster center and forms a binary 
with another star the cluster tends to becomes dominated by collisions. 
For an detailed description of this process we refer to Portegies Zwart et 
al. (1999). 

3.1. The collision rate from model calculations 

The number of collisions in each simulation ranges from 0-24, resulting a 
rate of 10~ 6 — 10~ 3 collision per star per Myr (see Tab.|2|)[]. These numbers 



1 Collision rates are computed over the time interval from zero age to the moment the 
last collision occurs. 
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are about two orders of magnitude higher than expected from simple cross 
section arguments. The reason for this discrepancy hides in the effect of 
mass segregation and the formation of binaries. 

The number of collisions depends on the number of stars and on the 
initial relaxation time of the stellar system. We therefore define f co ii as: 

Here ti ast is the time of the last collision. In Fig. [I] we give the collision 
rate f co u per star per million years as a function of the initial relaxation 
time of the model cluster. 

Models with the same initial conditions (Galactic tidal field, relaxation 
time and mass function) are averaged. 




Figure 1. Collision rate as function of the initial relaxation time 
for all models presented in Tab.|2|. The open circles give the results 
of our systems which are isolated from the Galactic potential (see 
Portegies Zwart et al 1999). Vertical bars represent Poissonian 
one-cr errors. The solid line is a least squares fit to the data. 

The solid line in Figgis a least squares fit to the results of our iV-body 
calculations, and yields 

K m n = 2.2 x 10- 4 t rbc -1 - [s -1 ]. (7) 
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The fit (Eq. |7|) is excellent, which is somewhat surprising since the 
models range a large area of parameter space, in initial density profile, 
mass function and tidal field (or not). Apparently these parameters hardly 
affect the collision rate. The fundamental parameters appear to be the 
initial relaxation time of the stellar system t T \ x and the number of stars N. 

3.2. The collision partners 

From cross section arguments one would expect that the most common 
collision candidates are among the stars which have the largest collective 
cross section (the collective projected area of the stars convolved with their 
gravitational focusing). Fig. Ogives the probability distribution of the stars 
which are most likely to be involved in a collision; the most likely collision 
counterparts are rather low mass stars, which are most common (see Porte- 
gies Zwart et al. 1999 for details). 




Figure 2. Distribution for collision candidates from derived 
using a Scalo (1986) mass function between 0.1 M and 100 M 
at zero age. 

Fig. |3] gives the probability distribution for the stars that did expe- 
rience a collision in our model calculations. The data are taken from a 
rather inhomogeneous sample of calculations (see Tab.0). 

The most striking difference between Fig.||| and Fig. [3] is the mean 
mass of the colliding stars. Where Fig. ^| only extends to 10 M stars and 
gives a highest probability for a collision between 0.7 M stars, Fig. |3| shows 
that much higher mass stars are much more likely to experience collisions. 
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The mean mass of the colliding stars according to the model calculations 
exceeds 10 M . The reasons for these discrepancies hide in the effect of 
mass segregation and the formation of binaries in the core of the stellar 
system (see Portegies Zwart et al. 1999 for details). 
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Figure 3. Probability distribution for collision candidates from 
our model calculations with Scalo (1986) initial mass function as 
in Fig.^. Darker shades indicates that more collisions between 
stars with these masses occurred in our calculations. 



3.3. Runaway growth 

The mean mass of a collision counterpart is about 10 M and many col- 
lisions involve the same star, the growth rate of this star is about 10 M 
per collision. Once the collision runaway explodes in a supernova or is 
ejected from the stellar system its growth stops and another star takes 
its place until the cluster dissolves in the tidal field of the Galaxy. (The 
models which were calculated without an external potential were stopped 
well before the cluster dissolved.) 

Fig. H presents the evolution tree for the collision sequence for two of 
the models. In this (an all other) cases a single object keeps colliding with 
other stars at a high rate. Only a few collisions occur between stars which 
do not finally coalesce with this collisions runaway. Only after the previous 
collision runaway is ejected from the cluster, another can take its place (see 
left panel in fig. |). 
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Figure 4. Coagulation tree for models R90W4 (left) and 
12k6A10 (right, see model Nk6X10 in Tab.|2|). Time is along 
the vertical axis from zero age (top) to the moment the last colli- 
sion occurs (bottom). The thickness of the lines scales linear with 
mass. 



4. Discussion 



In star clusters with a relaxation time smaller than a few 10 Myr collisions 
are quite common. The majority of collisions involve one selected object; 
the runaway merger. This seed object involves generally one of the initially 
most massive stars in the initial star cluster. If the cluster is embedded in 
the tidal field of the Galaxy, this selected star may eventually escape as a 
main sequence star or as the leftover from a supernova (a neutron star or 
black hole). Once the runaway merger has escaped another star will take 
its place until the cluster dissolves in the tidal field of the Galaxy. In an 
isolated cluster (or if the tidal fields on the cluster are small) it is harder 
to get rid of the collision runaway and in these cases the runaway merger 
may continue to grow in mass (see right panel in Fig.|3]). 

Most of our calculations, however, were computed with rather strong 
tidal forces, and these models dissolve within a few initial half mass relax- 
ation times. 

4.1. Observable characteristics of a collision runaway 

4.2. Super blue stragglers 

The stars which have experienced several collision may be brighter and 
bluer than any other star in the cluster, though detailed characteristics 
are hard to sketch. A few trivial characteristics may be those of a blue 
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straggler much more than twice the mass of the clusters' turn off (Sills 
& Lombardi 1997). For compact young star cluster such a characteristic, 
however, is not very useful as the turn off is ill defined for these systems. 

The three most massive stars in R136 have spectral type WN4.5 and 
appear to be younger than the other stars. Their age is estimated to be 
about lMyr (Massey & Hunter 1998; de Koter et al. 1997). These stars 
show violet absorption edges, which are common for late type (WN8 and 
later) stars but highly unusual for these early types (Conti et al. 1983). Also 
striking is that these stars are unusually hydrogen rich (Massey & Hunter 
1998). They are about an order of magnitude brighter than normal for such 
stars. Estimates for their masses range from 112 to 155 M (Chlebowski 
& Garmany 1991; Vacca et al. 1996). Two of them lie well inside the core 
of the cluster; the third is at a projected distance of about 0.6 pc from the 
core. 

Also the Quintuple cluster contains a very over-massive star, the Pistol 
star, with and estimated mass of about 150 M (Lang et al. 1999). Also 
this star appears to be rather odd and somewhat younger than the rest of 
the cluster stars. 

We aregue that both stars may be the result of (multiple) collisions, 
causing a supermassive star which is slightly younger due to the rejuvena- 
tion in the collision process. 

4.3. Super novae 

Mass segregation causes massive cluster members to be preferentially in the 
cluster center. This causes them to be excellent candidates for subsequent 
collisions. In our calculations the mean mass of a collision member is 
10 M Q . Stars with masses ^ 8M are expected to end their fuel burning 
life in a supernova, leaving a neutron star or a black hole. Collisions 
tend to reduce the number of type Ib/c and type II supernovae. Stars 
which normally would explode in a supernova (zero age mass ^ 8M ) are 
likely collision counterparts. A collision between two of such stars reduces 
therefore the number of supernova. 

4.4. Hypernovae? 

In compact young star clusters the mean accumulating mass (mean addi- 
tion in mass per collision) is about 10 M (see Fig. H). After A/" co ii collisions 
the initially most massive star (of 60 M to 100 M ) in the stellar system 
(the most likely collision runaway) has gained 10 Af co n M resulting in a 
total mass of 10A/" C oii + (60 — 100) M . Each collision may have rejuvenated 
the collision runaway somewhat, but eventually its core will grow until it 
exceeds a Chandrasekhar mass in iron, causing the star to explode in a 
supernova. The star itself is likely to be rapidly rotation, possibly close 
to break-up (Sills & Lombardi 1997). The enormous mass of the explod- 
ing star together with its high rotation rate may cause it to explode quite 
differently than stars in isolation. Nakamura et al. (2000) argues that 
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stars with an extremely high mass which are rotating very rapidly may 
be excellent candidates for Hypernovae, possibly even causing gamma-ray 
bursts. 

4.5. The formation of massive black holes 

As long as the collision runaway remains on the main sequence or one 
of the giant branches its cross section for further collisions remains large, 
due to the large size of the star. As soon as it explodes in a supernova 
things become somewhat uncertain as we do not understand the supernova 
mechanism very well. It is generally accepted that such massive stars form 
black holes but how much mass is lost from the collapsing star is uncertain. 
If the star loses a lot of mass in the supernova it may eject itself from the 
cluster, because the star was most likely to be a member of a close binary, 
which receives a high runaway velocity upon the supernova explosion (see 
Blaauw 1961; van den Heuvel et al. 2000; Portegies Zwart 2000). When 
little or no mass is lost in the supernova event the black hole remains in 
the cluster and may continue to grow in mass via subsequent collisions 
with other stars. Even though the size of the black hole has decreased 
dramatically it will still dominate the collision rate as its cross section is 
dominated by gravitational focusing. 

If we allowed little mass loss upon the supernova the formed black hole 
continues to grow in mass at about the same rate as before the supernova, 
indicating that the collision cross section is dominated by gravitational 
focusing. Take into account that if no mass is lost upon the supernova 
the black hole remains deep in the cluster core, where the collision are 
occurring. 

When more and more stars end their lifetime in a supernova the cluster 
will lose more mass. This results in a local expansion of the cluster core, 
finally terminating the collision efficiency We therefore do not expect that 
collisions remain dominating the cluster evolution for more than a few 10 
Myr. In older clusters binaries consisting of previously formed black holes 
start to heat the cluster (see Portegies Zwart & McMillan 2000). This 
will all finally result in a termination of the collision rate and causes the 
cluster to expand. Whether the cluster finally dissolves or remains bound to 
experience another phase of core collapse depends on details concerning the 
initial conditions of the star cluster. The further evolution of the cluster 
and its central massive black hole is a complicated matter and requires 
more study. We are quite happy that this matter is beyond the scope of 
this paper. 

4.6. The effect of primordial binaries 

All our calculations started without primordial binaries, which, of course, 
is an ill assumption. These clusters most likely contain a rich population 
of binaries, which have been formed together with the other stars in the 
cluster. For studying the collisional growth of a central object primordial 
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binaries might not be very important. Of course, the total mass in a binary 
is -on average- 1.5 times higher than the mean mass if single stars, but 
it will only make the increase in mass of the runaway merger larger and 
the absence of primordial binaries therefore reduces the effect of runaway 
growth. On the other hand, primordial binaries may heat the core at such a 
rate that collision are prevented altogether. This is very unlikely, as by the 
time binaries become hard enough to heat the cluster they are extremely 
vulnerable to collisions. Also, our runaway merger was generally a binary 
member, simply due to binary formation in the usual 3-body processes. 
We therefore expect the a rich population of primordial binaries does not 
lower the observed collision rate in our models. Detailed calculations in 
which we take a rich population of primordial binaries into account are in 
progress. 

5. Conclusions 

We performed a large number of A-body simulations which include the 
effects of stellar evolution, binary evolution and the tidal field of the Galaxy. 
The initial conditions are selected to represent the young and compact stars 
clusters such as Arches and the Quintuplet systems. 

Our model clusters are highly collisional, in the sense that collisions 
occur at a very high rate, much higher than expected from simple cross 
section arguments. The stars which are participating in collisions are gen- 
erally much more massive than the mean mass in the cluster. Generally 
one of the most massive stars in the cluster experiences multiple collisions 
until it either explodes in a supernova or is ejected from the cluster by a 
dynamical encounter. The collision rate per star per million years in our 
models scales very nicely with the initial relaxation time of the models, 
via: 

TZ coll = 2.2 x l(r 4 t rb r L0 [Myr- 1 ]. (8) 

Here t r \ x is the clusters' initial relaxation time in million years. 
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